function [L,NL,T] = trendClass(sigIn,SigSegmData,ParamWN)
%trendClass - Velocity trend classification before CNN training
%
%   [L,NL,T] = trendClass(sigIn,SigSegmData,ParamWN)
%
% The output L is a N-by-4 matrix L = [tRef,tIn,tFin,ClassTS], where
% tRef(k) is the reference time for the k-th velocity time series segment,
% tIn = tRef-Nc, tFin = tRef+Nd (see below) and ClassTS if the estimated
% class of this time series segment. 
% The output NL is the number of possible input levels (7 or 4, see below).
% The output T is a table whose variables are the possible classifications
% (including 'undefined' at the first row; therefore, the table has either
% 8 or 5 rows), the corresponding outcomes and the percentages. The table
% is also shown on the MATLAB Command Window.
%
% Input data:
%
%   sigIn      : input signal. It can be either a two columns matrix [t V]
%                or a three columns matrix [t R V], where t is the time,
%                R is the rainfall time series and  V is the velocity time
%                series. 
%
%   SigSegmData: cell variable provided by PluVelScalogram to generate the 
%                scalogram-bases images for the CNN training. This variable 
%                is 
%                   SigSegmData = {tRef,t1,t2,filena}.
%
% Moreover, these fields of ParamWN are used: 
%
%   Nc         : number of days to be considered before tRef; 
%
%   Nd         : number of days to be considered after tRef;
%
%   Vmh        : scalar/vector/matrix for the management of thresholds for 
%                the segments classification. Options:
%                - Vmh is a scalar. in this case, a 4-levels classification
%                  (including the transitios) is carried out, according to: 
%                       1: L, 2: L->H, 3: H, 4: H->L, 
%                  with the threshold VH = Vmh.
%                - Vmh is a 2-elements vector. In this case, a 7-levels
%                  classification is carried out with the thresholds 
%                  VM = Vmh(1) and VH = Vmh(2) (the vector is sorted in
%                  ascending order). The levels L, M and H are characterized
%                  by the conditions V<=M, M<V<=H and V>H respectively. The
%                  outputs are: 
%                       1: L, 2: L->M, 3: M, 4: M->H (or also L->H),
%                       5: H, 6: H->M, 7: M->L (or also H->L).
%                - Vmh is a 2-column matrix Vmh = [t VH], where VH(h) is
%                  the time-dependent threshold at t(h), as in the first
%                  case (four level output).
%                - Vmh is a 3-column matrix Vmh = [t VM VH], where VM(h) 
%                  and VH(h) are the time-dependent threshold at t(h), as 
%                  in the second case.
%
% If ParamWN is a char variable, it is supposed to be the name of a file
% having a field named 'ParamWN'. If the file loading or the extraction
% of the ParamWN value is not successful, it is ParamWN = [].
% If ParamWN is undefined or empty, a session of DefParam runs to generate
% ParamWN.
%
% See also PluVelInspect, PluVelScalogram, dataHomAug, dataTVT, 
%   TRLearnPluVel, DefParam.
%
%   [L,NL,T] = trendClass(sigIn,SigSegmData,ParamWN)

% G. Teza, 2020

if nargin < 3
    ParamWN = [];
end
if ischar(ParamWN)
    try 
        sv = load(ParamWN);
        ParamWN = sv.ParamWN;
    catch
        ParamWN = [];
    end
end
if isempty(ParamWN)
    ParamWN = DefParam;
end

if size(sigIn,2) == 2
    V = sigIn(:,2);
else 
    V = sigIn(:,3);
end
t = sigIn(:,1);

Nc = ParamWN.Nc;
Nd = ParamWN.Nd;

Vmh = ParamWN.Vmh;
if numel(Vmh) == 1
    Vt = [];
    Vm = [];
    Vh = Vmh;
    NL = 4;
elseif numel(Vmh) == 2
    Vt = [];
    Vmh = sort(Vmh);
    Vm = Vmh(1);
    Vh = Vmh(2);
    NL = 7;
elseif size(Vmh,2) == 2
    Vt = Vmh(:,1);
    Vm = [];
    Vh = Vmh(:,2);
    NL = 4;
elseif size(Vmh,2) == 3
    Vt = Vmh(:,1);
    Vm = Vmh(:,2);
    Vh = Vmh(:,3);
    NL = 7;    
else
    error('Invalid threshold input data');
end

nf = size(SigSegmData,1);

L = zeros(nf,4);
for k = 1:nf
    trefk = SigSegmData{k,1};
    L(k,1) = trefk;
    t1k = trefk-Nc;
    t2k = trefk+Nd;
    if (t1k >= t(1))&&(t2k <= t(end))
        L(k,2) = t1k;
        L(k,3) = t2k;
        [~,n1k] = min(abs(t-t1k));
        [~,n2k] = min(abs(t-t2k));
        if ~isempty(Vt)
            [~,nVt1k] = min(abs(Vt-t1k));
            [~,nVt2k] = min(abs(Vt-t2k));
            Vhk = mean(Vh(nVt1k:nVt2k));
            if ~isempty(Vm)
                Vmk = mean(Vm(nVt1k:nVt2k));
            else
                Vmk = [];
            end
        else
            Vmk = Vm;
            Vhk = Vh;
        end
        Vk = V(n1k:n2k);
        lhk = round(0.5*length(Vk));
        Vk1 = Vk(1:lhk);
        Vk2 = Vk(lhk:end);
        VNk1 = higherN(Vk1,5);
        VNk2 = higherN(Vk2,5);
        if ~isempty(Vmk)
            % CASE 7 OUTPUTS
            if (nanmean(VNk1)<=Vmk) && (nanmean(VNk2)<=Vmk) 
                L(k,4) = 1;     % L -> L
            elseif (nanmean(VNk1)<=Vmk) && (nanmean(VNk2)>Vmk) && (nanmean(VNk2)<=Vhk)
                L(k,4) = 2;     % L -> M
            elseif (nanmean(VNk1)<=Vmk) && (nanmean(VNk2)>Vhk)  % direct 1-5 transition
                L(k,4) = 4;     % L -> H
            elseif (nanmean(VNk1)>Vmk) && (nanmean(VNk1)<=Vhk) && (nanmean(VNk2)>Vmk) && (nanmean(VNk2)<=Vhk)
                L(k,4) = 3;     % M -> M
            elseif (nanmean(VNk1)>Vmk) && (nanmean(VNk1)<=Vhk) && (nanmean(VNk2)>Vhk)
                L(k,4) = 4;     % M -> H
            elseif (nanmean(VNk1)>Vhk) && (nanmean(VNk2)>Vhk)
                L(k,4) = 5;     % H -> H
            elseif (nanmean(VNk1)>Vhk) && (nanmean(VNk2)>Vmk) && (nanmean(VNk2)<=Vhk)
                L(k,4) = 6;     % H -> M
            elseif (nanmean(VNk1)>Vmk) && (nanmean(VNk1)<=Vhk) && (nanmean(VNk2)<=Vmk)
                L(k,4) = 7;     % M -> L
            elseif (nanmean(VNk1)>Vhk) && (nanmean(VNk2)<=Vmk)
                L(k,4) = 6;     % H -> L    
            end
        else
            % CASE 4 OUTPUTS
            if (nanmean(VNk1)<=Vhk) && (namean(VNk2)<=Vhk) 
                L(k,4) = 1;     % L -> L
            elseif (namean(VNk1)<=Vhk) && (namean(VNk2)>Vhk)
                L(k,4) = 2;     % L -> H
            elseif (namean(VNk1)>Vmk) && (namean(VNk2)>Vhk)  
                L(k,4) = 3;     % H -> H
            elseif (namean(VNk1)>Vhk) && (namean(VNk1)<=Vhk)
                L(k,4) = 4;     % H -> L
            end
        end
    end
end

if NL == 7
    CL = {'Undefined','L1','L2','L3','L4','L5','L6','L7'}';
    VL = zeros(8,1);
    for k = 1:8
        Ik = L(:,4) == k-1;
        VL(k) = sum(Ik);
    end
 else
    CL = {'Undefined','L1','L2','L3','L4'}';
    VL = zeros(5,1);
    for k = 1:5
        Ik = L(:,4) == k-1;
        VL(k) = sum(Ik);
    end
end
PL = VL/nf*100;

T = table(CL,VL,PL,'VariableNames',{'classification','number','percentage'});
fprintf('\nSummary about time series classification:\n\n'); 
disp(T);


%% Ancillary functions
function VN = higherN(V,N)
% VN is the set of the N higher elements of V
VN = zeros(1,N);
Vk = V;
for k = 1:N
    [maxk,Ik] = max(Vk);
    VN(k) = maxk;
    Vk(Ik) = [];
end       